--- layout: page title: PO2 gradiants --- diffusion_advection_regressions
In [19]:
import numpy as np
import pandas as pd

import statsmodels.api
import scipy.stats

import bokeh.io
import bokeh.plotting

#local .py file for some plotting functions and non-parametric bootstrapping utils
import plotting_utils

import numba

bokeh.io.output_notebook()
Loading BokehJS ...

We can now read in the data into a dataframe for analyis.

In [23]:
df = pd.read_csv("./20190322_supp_table_2.csv")

We will add some calculations/generate a species averaged version of the dataframe.

In [25]:
df['species_underscore'] = [spec.replace(" ", "_") for  spec in df['species']]
df_averages = df.groupby(['species', 'species_underscore', 'spiracle'], as_index=False).aggregate(np.average)
df_averages['subfamily'] = df.groupby(['species', 'species_underscore', 'spiracle'], as_index=False).aggregate(max)['subfamily']
species_per_subfam=df_averages.groupby(['subfamily', 'spiracle'], as_index=False).count().groupby('subfamily').aggregate(max).reset_index()[['subfamily', 'species']]
species_per_subfam.columns = ('subfamily', 'subfam_count')
df_averages = df_averages.merge(species_per_subfam, on='subfamily')

df_averages['log area (mm^2)'] = np.log10(df_averages['area (mm^2)'])
df_averages['log dist'] = np.log10(df_averages['depth (mm)'])
df_averages['log mass (g)'] = np.log10(df_averages['mass (g)'])
df_averages['log area/dist'] = np.log10(df_averages['area (mm^2)']/df_averages['depth (mm)'])
df_averages['log area^2/dist'] = np.log10(df_averages['area (mm^2)']**2/df_averages['depth (mm)'])

df['log area (mm^2)'] = np.log10(df['area (mm^2)'])
df['log dist'] = np.log10(df['depth (mm)'])
df['log mass (g)'] = np.log10(df['mass (g)'])
df['log area/dist'] = np.log10(df['area (mm^2)']/df['depth (mm)'])
df['log area^2/dist'] = np.log10(df['area (mm^2)']**2/df['depth (mm)'])

df_averages.head()
Out[25]:
species species_underscore spiracle Unnamed: 0 area (mm^2) depth (mm) log area (mm^2) log area/dist log area^2/dist log dist log mass (g) mass (g) subfamily subfam_count
0 Coelorrhina hornimani Coelorrhina_hornimani 1 87.0 0.135347 0.416717 -0.868551 -0.488392 -1.356943 -0.380159 0.053078 1.13 Cetoniinae 6
1 Coelorrhina hornimani Coelorrhina_hornimani 2 70.0 0.084207 0.451409 -1.074651 -0.729221 -1.803872 -0.345430 0.053078 1.13 Cetoniinae 6
2 Coelorrhina hornimani Coelorrhina_hornimani 3 53.0 0.106693 0.325444 -0.971862 -0.484339 -1.456201 -0.487524 0.053078 1.13 Cetoniinae 6
3 Coelorrhina hornimani Coelorrhina_hornimani 4 36.0 0.115574 0.481558 -0.937142 -0.619790 -1.556932 -0.317351 0.053078 1.13 Cetoniinae 6
4 Coelorrhina hornimani Coelorrhina_hornimani 5 19.0 0.119145 0.506751 -0.923923 -0.628717 -1.552640 -0.295205 0.053078 1.13 Cetoniinae 6

Using these morphological properties of the spiracles that we measured, we can calculate real world meaningful gas exchange properties with some unit conversions/physical gas properties. To calculate the diffusive capacity of the spiracles, we use the following relation:

\begin{align} G_{\mathrm{diff}} = \frac{\text{area}}{\text{depth}} * D * \beta \end{align}

with $D$ (the diffusivity constant for O2 in air) $= 0.178 \text{ cm}^2$ and $\beta$ (the capacitance coefficient for oxygen in air) $= 404 \frac{\text{ nmol}}{\text{cm}^{3} \text{kPa}}$. With that value, you can calculate the oxygen consumption rate in $\frac{\text{nmol}}{\text{sec}}$ by multiplying the oxygen partial pressure gradient $\Delta \text{PO}_2$ across the spiracle in kPa by $G_{\text{diff}}$, that is

\begin{align} \text{Oxygen consumption rate }\left(\frac{\text{nmol}}{\text{sec}}\right) = \Delta \text{PO}_2 * G_{\mathrm{diff}} \end{align}

To calculate advective gas transport capacities, we use $G_{\text{adv}}$, which comes from Poiseulle’s law and is

\begin{align} G_{\mathrm{adv}} = \frac{\text{area}^2}{8 * \pi * \text{dynamic viscosity} * \text{depth}} \end{align}

where the dynamic viscosity of air is $1.86 * 10^{-8} \text{kPa sec}$. With these relations in hand, we can calculate both $G_{\text{diff}}$ and $G_{\text{adv}}$ for the beetle spiracles, both individually, and also by summing the individual capacities for all spiracles in the beetle.

In [9]:
df_test = df_averages.copy()
df_test['area/depth'] = (df_averages['area (mm^2)']/df_averages['depth (mm)'])/10
df_test['area^2/depth'] = ((df_averages['area (mm^2)']**2/df_averages['depth (mm)']))/(1000*1000*1000) #convert to cubic meters

df_test['g_diff'] = df_test['area/depth']*0.178*404
df_test['g_adv'] = df_test['area^2/depth']*(1/(np.pi*8*1.86*10**(-8)))

df_summed = df_test.groupby('species').median().reset_index()[['species', 'log mass (g)', 'g_diff', 'g_adv']]
df_summed['g_diff'] = df_test.groupby('species').sum().reset_index()['g_diff']
df_summed['g_adv'] = df_test.groupby('species').sum().reset_index()['g_adv']
In [27]:
plots = []
resample_size = 10_000
lw = 2
cs = 12

m = df_summed['log mass (g)']
g_diff = np.log10(df_summed['g_diff']*2) #double the number for g_diff since we only measured one side of the animal, they have two of each spiracle

slope, intercept = np.polyfit(m, g_diff, deg=1)
x = np.array([m.min(), m.max()])
y = slope * x + intercept

p = bokeh.plotting.figure(plot_height=300, plot_width=400,
                          x_range=(x[0]-0.1, x[1]+0.1), y_range=(g_diff.min()-0.1, g_diff.max()+0.2))
p.outline_line_color = None
p.yaxis.minor_tick_line_color = None
p.xaxis.minor_tick_line_color = None
p.grid.grid_line_color = None

slope_comp = 0.75
intercept1 = plotting_utils.first_intercept(slope_comp, x.max(), g_diff.min()) -0.5
line_scale = (y.max() - y.min())/5
around_line=0.2
for i in line_scale*np.array(range(30))+intercept1:
        try:
            lx, ly = plotting_utils.generate_line(intercept=i, slope=slope_comp,
                                                   bounds=(x[0]-around_line, x[1]+around_line,
                                                           g_diff.min()-around_line, g_diff.max()+around_line+0.4), point=x[1])
            p.line(lx, ly, color='grey', alpha=0.3)
        except:
            pass
        
bs_slope_reps, bs_intercept_reps, _, _ = plotting_utils.draw_bs_pairs_linreg(m, g_diff, size=resample_size)
p.title.text = ('Slope: '+ str(round(slope, 2)) + ' intercept: ' + str(round(intercept, 2)) +' slope 95% CI: ' + 
                str([round(j, 3) for j in np.percentile(bs_slope_reps, [2.5, 97.5])]) + ' ' + str(np.sum(bs_slope_reps > 0.75)))
x_boot = np.linspace(m.min(), m.max(), 200)
y_boot = np.outer(bs_slope_reps, x_boot) + np.stack([bs_intercept_reps]*200, axis=1)
low, high = np.percentile(y_boot, [2.5, 97.5], axis=0)
p1 = np.append(x_boot, x_boot[::-1])
p2 = np.append(low, high[::-1])
p.patch(p1, p2, alpha=0.5, color='lightgrey')


p.circle(m, g_diff, color='black', size=cs)
p.line(x, y, color='black', line_width=lw, line_cap='round')
p.output_backend='svg'
plots.append(p)
#bokeh.io.show(p)

m = df_summed['log mass (g)']
g_adv = np.log10(df_summed['g_adv']*2) #double the number for g_diff since we only measured one side of the animal, they have two of each spiracle

slope, intercept = np.polyfit(m, g_adv, deg=1)
x = np.array([m.min(), m.max()])
y = slope * x + intercept

p = bokeh.plotting.figure(plot_height=300, plot_width=400, x_range=(x[0]-0.1, x[1]+0.1), y_range=(g_adv.min()-0.2, g_adv.max()+0.5))
p.outline_line_color = None
p.yaxis.minor_tick_line_color = None
p.xaxis.minor_tick_line_color = None
p.grid.grid_line_color = None

slope_comp = 0.75
intercept1 = plotting_utils.first_intercept(slope_comp, x.max(), g_adv.min())
line_scale = (y.max() - y.min())/7
around_line=0.4
for i in line_scale*np.array(range(30))+intercept1:
        try:
            lx, ly = plotting_utils.generate_line(intercept=i, slope=slope_comp,
                                                   bounds=(x[0]-around_line, x[1]+around_line,
                                                           g_adv.min()-around_line, g_adv.max()+around_line+0.5), point=x[1])
            p.line(lx, ly, color='grey', alpha=0.3)
        except:
            pass


bs_slope_reps, bs_intercept_reps, _, _ = plotting_utils.draw_bs_pairs_linreg(m, g_adv, size=resample_size)
p.title.text = ('Slope: '+ str(round(slope, 2)) + ' intercept: ' + str(round(intercept, 2)) +' slope 95% CI: ' + 
                str([round(j, 3) for j in np.percentile(bs_slope_reps, [2.5, 97.5])]) + ' ' + str(np.sum(bs_slope_reps < 0.75)))
x_boot = np.linspace(m.min(), m.max(), 200)
y_boot = np.outer(bs_slope_reps, x_boot) + np.stack([bs_intercept_reps]*200, axis=1)
low, high = np.percentile(y_boot, [2.5, 97.5], axis=0)
p1 = np.append(x_boot, x_boot[::-1])
p2 = np.append(low, high[::-1])
p.patch(p1, p2, alpha=0.5, color='lightgrey')

p.circle(m, g_adv, color='black', size=cs)
p.line(x, y, color='black', line_width=lw, line_cap='round')
p.output_backend='svg'
plots.append(p)

bokeh.io.show(plots[0])
bokeh.io.show(plots[1])

To calculate the ∆PO2 across the spiracles needed to supply beetle metabolic demand by diffusion, the metabolic rate for a quiescent beetle at a body temperature of 25 °C of a given mass was estimated from here with the following equation: $\text{log}_{10} (\text{metabolic rate } \mu \text{W}) = 3.20 + 0.75*\mathrm{log}_{10}(\text{mass (g)}$ and the assumption of 20.7 kJ/L for O2. For O2 at 25 °C, 24.5 mol/L was also assumed. Based on here, flight metabolic rate of small insects is in the range of 8x resting, whereas it is on the order of 32x resting metabolic rate in large insects. Most scarab beetles are endothermic during flight, so flight metabolic rates of these warm beetles could double this value to 64x with a 10 °C increase in thoracic temperature if the Q10 is 2. Even this calculation may be a conservative estimate of hovering metabolic rate, as oxygen consumption rose to 90x higher than those of quiescent, 25 °C fig beetles in one study, and maximal flight metabolic rates may be 1.25-2x higher than during hovering. Therefore, we estimated maximal aerobic metabolic rate during flight as 90x those of quiescent beetles. The required PO2 gradient across the spiracles to support gas exchange by diffusion at rest and during flight was calculated by rearranging our previous $G_{text{diff}}$ equation and performing unit conversions as follows:

\begin{align} \Delta\text{PO}_2 (\text{kPa}) = \left( \frac{\left(10^{3.20 + 0.75*\mathrm{log}_{10}(\text{mass (g)})}\mu W\right) \left( \frac{1 \frac{\mu J}{s}}{1 \mu W} \right)\left( \frac{nL}{20.7 \mu J} \right) \left( \frac{nmol}{24.5 nL} \right)}{\left(\frac{\text{area}}{\text{depth}}(cm)\right) \left( 0.178 \frac{cm^2}{sec}\right) \left( 404 \frac{nmol}{cm^3 kPa}\right)}\right) \end{align}

With this relation in hand, we plot the required oxygen partial pressure gradient across just the spiracle needed to supply the metabolic demand of the insect.

In [28]:
m = df_summed['log mass (g)']
g_diff = np.log10(df_summed['g_diff']*2)
slope, intercept = np.polyfit(m, g_diff, deg=1)
po2 = np.log10((1*((10**(3.20 + 0.75*m))*(1/20.7)*(1/24.5))/(df_summed['g_diff']*2)))
print(slope)
N = 200
y_top = 90
x = np.linspace(m.min(), m.max(), N)
y = np.linspace(1, y_top, N)
im = np.zeros((N, N))
for j, yi in zip(range(N), y):
    for i, xi in zip(range(N), x):
        im[i, j] = (yi*((10**(3.20 + 0.75*xi))*(1/20.7)*(1/24.5))/(10**(slope * xi + intercept)))

        
N = 800
y_top = 90
x = np.linspace(m.min(), m.max(), N)
y = np.linspace(1, y_top, N)

x_pos = []
for j, yi in zip(range(N), y):
    [x_pos.append([xj, yj, zj]) for xj, yj, zj in zip(x, np.log10(yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/(10**(slope * x + intercept))),
                                                                  yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/(10**(slope * x + intercept)))]
    
x_min, y_min, v_min = (np.array(x_pos)[:, 0].min(), np.array(x_pos)[:, 1].min(), np.array(x_pos)[:, 2].min())
x_stride, y_stride = ((np.array(x_pos)[:, 0].max() - x_min)/N, (np.array(x_pos)[:, 1].max() - y_min)/N)

im = np.ones((N, N))*v_min
for xj, yj, vj in x_pos:
    im[int(np.ceil((yj-y_min)/y_stride-1)), int(np.ceil((xj-x_min)/x_stride-1))] = vj
im = scipy.ndimage.gaussian_filter(im, sigma=1)
            
p = bokeh.plotting.figure(tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")], plot_height=300, plot_width=400)
#p.x_range.range_padding = p.y_range.range_padding = 0

cmap = bokeh.models.LinearColorMapper(palette='Viridis256', low=im.min(), high=im.max())
cmap_low, cmap_high = (np.min(1*(10**po2)), np.max(90*(10**po2)))
cmap = bokeh.models.LinearColorMapper(palette='Viridis256', low=cmap_low, high=cmap_high)

p.image(image=[im], x=x.min(), y=np.log10(im.min()), dw=x.max()-x.min(), dh=np.log10(im.max())-np.log10(im.min()), color_mapper=cmap, level="image", )

color_bar = bokeh.models.ColorBar(color_mapper=cmap, location=(0,0), ticker=bokeh.models.BasicTicker(desired_num_ticks=12, base=10))
p.add_layout(color_bar, 'right')
p.grid.grid_line_color = None

p.patch([x.min(), x.max(), x.max(), x.min()], [np.log10((8*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).min(),
                                               np.log10((8*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).max(),
                                               np.log10(im.min()), np.log10(im.min())], color='white', line_width=2)

p.patch([x.min(), x.max(), x.max(), x.min()], [np.log10((90*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).min(),
                                               np.log10((90*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).max(),
                                               np.log10(im.max()), np.log10(im.max())], color='white', line_width=2)

p.line([x.min(), x.max()], [np.log10(21), np.log10(21)], color='lightgrey', line_width=2, alpha=0.75)

#p.patch([x.min(), x.max(), x.max(), x.min()], [np.log10((1*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).min(),
#                                               np.log10((1*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).max(),
#                                               np.log10((90*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).max(),
#                                               np.log10((90*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))).min()], color='white', line_width=0, alpha=0.2)


line_color='black'
#[p.line(x, np.log10((yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))), color=line_color) for yi in np.linspace(1, y_top, 10)]
[p.line(x, np.log10((yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(slope * x + intercept))), color=line_color, line_width=lw, line_cap='round') for yi in [1, 8, 90]]
#[p.line(x, np.log10((yi*((10**(3.20 + 0.75*x))*(1/20.7)*(1/24.5))/10**(0.33  * x + intercept))), color='white', line_dash='dashed', line_width=5, line_cap='round') for yi in [1, 8, 90]]

dot_color='white'
dot_line='black'
dot_width = 1
cmapper = bokeh.transform.linear_cmap('c', palette='Viridis256', low=np.min(1*(10**po2)), high=np.max(90*(10**po2)))
for mult in [1, 8, 90]:
    source = bokeh.models.ColumnDataSource(data=dict(x=m, y=np.log10(mult*(10**po2)), c=mult*(10**po2)))
    p.circle('x', 'y', source = source, size=cs, line_color=dot_line, line_width = dot_width, fill_color=cmapper)
p.outline_line_color = None
p.yaxis.minor_tick_line_color = None
p.xaxis.minor_tick_line_color = None
p.output_backend='svg'
plots.append(p)

bokeh.io.show(bokeh.layouts.gridplot([plots[0], plots[2], plots[1]], ncols=3))

0.39370488659751474

The required partial pressure gradient with a 90x resting metabolic demand is

In [31]:
90*(10**po2)
Out[31]:
0     8.047311
1     4.574343
2    17.376828
3    41.177161
4    13.268803
5    30.720956
6    21.797164
7     9.780347
8    10.284551
9    41.391985
dtype: float64

These values exceed the partial pressure of air in the atmosphere! No way a large beetle can use diffusion for gas exchange, it must use active advective processes.

In [29]:
#TO EXPORT THE PLOTS FROM ABOVE
#bokeh.io.export_svgs(plots[2], filename="./ptest.svg")
In [15]:
%reload_ext watermark
%watermark -p bokeh

bokeh 1.4.0
In [ ]: